Content under Creative Commons Attribution license CC-BY 4.0, code under MIT license (c)2014 P. Y. Chuang.

Potential Flow Over an Airfoil with Finite Difference Method

Introduction

Potential-flow theory has been used to study the flow over 2D airfoils for a long time. Using complex analysis, we can obtain analytical solutions for many potential flows, especially the flow over various 2D airfoils. When computers became available, people could also solve potential flows numerically using panel methods, finite-difference methods, boundary-element methods, etc.

Though there are many unrealistic assumptions in potential flow, the theory provides a good beginning for further analyses. For example, some CFD codes can solve potential flow first, in order to provide a good initial guess of viscous flow solver.

In 2D potential flow, we can solve for either the potential function or stream function through their governing equtaions: the Laplace equation. The governing equation for stream function is $$ \nabla^2\Psi=\frac{\partial^2 \Psi}{\partial x^2}+\frac{\partial^2 \Psi}{\partial y^2}=0 $$ where $\Psi$ is the stream function. If the velocity vector $\vec{V}=u\vec{i}+v\vec{j}$ represents the flow velocity, the relationship between flow velocity and the stream function is: $$ u=\frac{\partial \Psi}{\partial y}\\ v=-\frac{\partial \Psi}{\partial x} $$

In this notebook, we show how to solve potential flow over 2D airfoils using the finite difference method.

NACA 4-digit Airfoil Series

We use NACA 4-digit airfoils for our example in this notebook. For an introduction to NACA 4-digit airfoils, refer to reference [1].

Let's first define a Python function that can be used to generate the thickness profile of a NACA airfoil, given by：

$$ y_t = 5t\left[ 0.2969\sqrt{\frac{x}{c}} -0.126\left(\frac{x}{c}\right) -0.3516\left(\frac{x}{c}\right)^2 +0.2843\left(\frac{x}{c}\right)^3 -0.1036\left(\frac{x}{c}\right)^4 \right] $$

where $y_t$, $t$, $x$, and $c$ represent the thickness, maximum thickness, $x$-coordinate and chord length, respectively. （Note: the last coefficient is modified from $-0.1015$ to $-0.1036$ in order to "close" the trailing edge. If we used $-0.1015$, there would be a thickness at the trailing edge, just like airfoils in the real world. But in the computational world, we prefer to use zero-thickness trailing edges.)

Note that $y_t$ is in fact just the half thickness of the airfoil, so the upper and lower profiles must be generated separately. The lower part takes the negetive values of those corresponding to the upper part. Using this function, we can generate the profiles for airfoils with the first and the second digits being zero (i.e., symmetrical airfoils). Here we show the profile of the NACA-0014 airfoil as an example:

For cambered 4-digit arifoils (i.e., the first and the second digit are not zero), we need to define Python functions for center lines (mean camber lines) and the angles, $\theta$, between center lines and the horizontal line.

$$ y_c = \left\{\begin{array}{lr} m\frac{x}{p^2}\left(2p-\frac{x}{c}\right)\text{,} && 0\le x\le p\times c\\ m\frac{c-x}{(1-p)^2}\left(1+\frac{x}{c}-2p\right)\text{,} && p\times c \le x \le c \end{array}\right. $$

and

$$ \theta=\tan^{-1}{\frac{dy_c}{dx}} $$

where

$$ \frac{dy_c}{dx}=\left\{\begin{array}{lr} \frac{2m}{p^2}\left(p-\frac{x}{c}\right)\text{,} && 0\le x\le p\times c\\ \frac{2m}{(1-p)^2}\left(p-\frac{x}{c}\right)\text{,} && p\times c \le x \le c \end{array}\right. $$

Finally, the relationship between the coordinates on the upper and lower surfaces of cambered airfoils ($x_U$, $y_U$, $x_L$, and $y_L$), the center line coordinate ($x$ and $y_c$) and the thickness ($y_t$) are $$ \begin{array}{ll} x_U=x-y_t\sin{\theta}\text{,} && y_U=y_c+y_t\cos{\theta}\\ x_L=x+y_t\sin{\theta}\text{,} && y_L=y_c-y_t\cos{\theta} \end{array} $$

We wrote a function that can generate both symmetrical and cambered NACA 4-digit airfoils and use a sign to represent upper or lower profile:

The following shows the profile of a NACA 2412 airfoil as an example.

Discretization in Space of NACA 2412 Airfoil

Physical Domain and Computational Doamin

Here we introduce the concepts of physical and computational domains, following chapter 9 of Hoffmann & Chiang, 2000 [2].

With the finite difference method, rectangular grids work best. But it is difficult to discretize the domain into rectangular grids when there are objects (like airfoils) within it.

In such a case, the domain can be discretized into a non-Cartesian structured mesh. Such a mesh may contain arbitrary quadrilateral grids, instead of rectangular grids. The approach is to map our real domain into another space that is discretized into rectangular grids and solve our problem in this space. The original domain is called physical domain, and the space in which we solve our problems is called computational domain. The following figure shows the concepts of physical and computational domains.

Image credit: Hoffmann & Chiang, 2000 (Fig. 9-23a and b).

We define the coordinates in the computational domain to be $\xi$ and $\eta$ (called computational coordinates) and those in the physical doamin to be $x$ and $y$ (physical coordinates). Physical coordinates are functions of computational coordinates, and vice versa. This means that all calculations performed in the physical domain can also be performed in the computational domain through the relationship between $x$, $y$ and $\xi$, $\eta$, after their relationship has been determined.

If we fix the computational coordinates of grid points in a way that all grids in computational domain are rectangular grids with $\Delta\xi=\Delta\eta=1$, what we need to do next is to determine the physical coordinates of these computational grid points. Let $i$ and $j$ be the indices of grid points in the computational domain. What we want are the values $x_{i,j}, y_{i,j}$ representing the $x, y$ components in physical space of the grid points $(i,j)$ in the computational domain.

In this notebook, we use an O-type mesh on the physical domain, that is, the outer boundary of the physical domain is a circle, as shown in the figure below.

Image credit: Hoffmann & Chiang, 2000 (Fig. 9-28).

We divide the annular area between the arifoil and the outer boundary by a line $\overline{AC}$:

Image credit: Hoffmann & Chiang, 2000 (Fig. 9-29).

Then we "expand" the divided annular area to a rectangular domain, which will be our computational domain:

Image credit: Hoffmann & Chiang, 2000 (Fig. 9-30).

From the name of boundaries ($B_1$~$B_4$), we can understand the relationship between our physical domain and the computational domain.

Grid Generation

The grids on our computational domain are simply rectangular grids. The lengths of boundaries and the size of the grids in the computational domain are not important, because they only make a difference on the transformation relationships between two domains. So let's fix the size of grids in our computational domain to 1, i.e., $\Delta\xi=\Delta\eta=1$. Also, the numbers of nodes on the $\xi$ and $\eta$ directions are $N_{\xi}=51$ and $N_{\eta}=21$. The grids on the computational domain are defined in code as follows:

The next step is to determine the physical coordinates of the grid points, $x_{i,j}$ and $y_{i,j}$, in the computational domain. In the following two subsections, we introduce two grid generation methods: the algebraic method and elliptic method.

Algebraic Grid Generation

The idea of algebraic grid generation is to interpolate the physical coordinates of interior grid points between two known boundaries, i.e., the airfoil and the outer boundary in our case.

We know that the number of grid points on the airfoil and outer boundary should be the same with $N_{\xi}$, while on the dividing line, $\overline{AC}$, they should equal $N_{\eta}$. We also know the physical coordinates of grid points on $\eta=0$ and $\eta=N_{\eta}$, because they are located on the airfoil surface and the outer boundary and can be determined as follows:

$$ \begin{array}{ll} x_{i, j=0} = x_{airfoil}\text{,} && y_{i, j=0} = y_{airfoil}\\ x_{i, j=-1} = x_{outer\ boundary}\text{,} && y_{i, j=-1} = y_{outer\ boundary} \end{array} $$

For convenience, we use Python indexing, where $j=-1$ represents the last grid points in the $\eta$ direction.

The $x$ and $y$ coordinates for interior points (i.e., $0\le i\le -1$ and $1\le j \le -2$, using Python indexing) can be obtained through interpolation. This is called algebraic grid-generation.

We now generate algebraic grids for a NACA 2412 airfoil. The ordering of grid points on the airfoil is clockwise and starts from the trailing edge. The center of the circular outer boundary is located at $0.5c$.

Now, we interpolate the coordinates of interior grid points.

Let's see what the mesh looks like.

Looking closer to the vicinity of the airfoil:

Though the mesh quality wasn't actually examined here, it does not looks good near the airfoil. That's why we need another grid-generation method.

Elliptic Grid Generation

In the elliptic grid-generation method, we use elliptic PDEs to constrain the relationship between physical and computational coordinates. The elliptic PDEs used in this grid generation method are:

$$ \frac{\partial^2 \xi}{\partial x^2}+\frac{\partial^2 \xi}{\partial y^2}=0\\ $$$$ \frac{\partial^2 \eta}{\partial x^2}+\frac{\partial^2 \eta}{\partial y^2}=0\\ $$

However, if we solve these PDEs (if we could), what we would obtain are the computational coordinates, i.e., $\xi=\xi(x, y)$ and $\eta=\eta(x, y)$, which are not what we want, because we already fixed the computational coordinates (they're simply on a rectangular grid with $\Delta\xi=\Delta\eta=1$).

Therefore, the actual PDEs that we are going to solve, which are derived from the above PDEs (see reference [2]), are

$$ a\frac{\partial^2 x}{\partial \xi^2}-2b\frac{\partial^2 x}{\partial \xi \partial \eta}+c\frac{\partial^2 x}{\partial \eta^2}=0 $$$$ a\frac{\partial^2 y}{\partial \xi^2}-2b\frac{\partial^2 y}{\partial \xi \partial \eta}+c\frac{\partial^2 y}{\partial \eta^2}=0 $$

where

$$ a=\left(\frac{\partial x}{\partial \eta}\right)^2+\left(\frac{\partial y}{\partial \eta}\right)^2\\ b=\left(\frac{\partial x}{\partial \xi}\right)\left(\frac{\partial x}{\partial \eta}\right)+\left(\frac{\partial y}{\partial \xi}\right)\left(\frac{\partial y}{\partial \eta}\right)\\ c=\left(\frac{\partial x}{\partial \xi}\right)^2+\left(\frac{\partial y}{\partial \eta}\right)^2 $$

We can see that, in order to solve the PDEs, we need two boundary conditions on each direction, This means that we need $x_{i,j=0}$, $x_{i,j=-1}$, $x_{i=0,j}$, $x_{i=-1,j}$, $y_{i,j=0}$, $y_{i,j=-1}$, $y_{i=0,j}$, and $y_{i=-1,j}$.

Luckly, $x_{i,j=0}$, $x_{i,j=-1}$, $y_{i,j=0}$, and $y_{i,j=-1}$ are known and are the physical coordinates of nodes on the airfoils and outer boundary, as described in the previous subsection.

From the illustration above, we also know that $x_{i=0,j}$ and $x_{i=-1,j}$ are actually the same, because they both represent the physical coordinates of nodes on the dividing line $\overline{AC}$. Similarly, $y_{i=0,j}$ and $y_{i=-1,j}$ are the same, too. This implies that periodic boundary conditions can be applied on $x_{i=0,j}$, $x_{i=-1,j}$, $y_{i=0,j}$, and $y_{i=-1,j}$.

These PDEs are non-linear but we linearize the solution by letting the coefficients lag behind by one iteration. We'll use the basic Jacobi iterative method with the results obtained from the algebraic grid-generation method as initial guess.

We can use central differences on these PDEs, and with $\Delta\xi=\Delta\eta=1$, the discretized equations are very simple. Since the forms of these two PDEs are the same, their discretized PDEs are also the same. The discretized equation of $x$ are

$$ a(x_{i+1, j}-2x_{i,j}+x_{i-1, j}) - \frac{1}{2}b(x_{i+1, j+1}-x_{i+1, j-1}+x_{i-1,j-1}-x_{i-1,j+1})+c(x_{i, j+1}-2x_{i,j}+x_{i, j-1})=0 $$

or

$$ x_{i,j}=\frac{1}{2}\left\{ a(x_{i+1, j}+x_{i-1, j}) - \frac{1}{2}b(x_{i+1, j+1}-x_{i+1, j-1}+x_{i-1,j-1}-x_{i-1,j+1})+c(x_{i, j+1}+x_{i, j-1}) \right\} / (a+c) $$

where

$$ a=(\frac{x_{i,j+1}-x_{i,j-1}}{2})^2+(\frac{y_{i,j+1}-y_{i,j-1}}{2})^2\\ b=(\frac{x_{i+1,j}-x_{i-1,j}}{2})(\frac{x_{i,j+1}-x_{i,j-1}}{2})+(\frac{y_{i+1,j}-y_{i-1,j}}{2})(\frac{y_{i,j+1}-y_{i,j-1}}{2})\\ c=(\frac{x_{i+1,j}-x_{i-1,j}}{2})^2+(\frac{y_{i+1,j}-y_{i-1,j}}{2})^2 $$

Replace $x$ with $y$ when we want to solve $y$.

We can now solve these equations. The criterion for stopping the iteration is that the maximum difference between a current result and the result of the previous iteration should be less than $10^{-6}$.

Let's see what the mesh looks like.

The grids near the airfoil look more clustered. This is good. And how about the grids in the vicinity of the airfoil?

It looks better than what it looks like in algebraic grid generation method. We can use it to solvo our potential flow now!

Potential Flow over NACA 2412

Governing Equation in the Computational Domain

We know the governing equation of the stream function $\Psi$ in the physical domain is the Laplace equation, as described in the introduction section. To solve the stream function in the computational domain, we need a new governing equation. According to the derivation in reference [3], we obtain the governing equation for stream function in the computational domain as follows:

$$ a\frac{\partial^2 \Psi}{\partial \xi^2}-2b\frac{\partial^2 \Psi}{\partial \xi \partial \eta}+c\frac{\partial^2 \Psi}{\partial \eta^2}=0 $$

where

$$ a=\left(\frac{\partial x}{\partial \eta}\right)^2+\left(\frac{\partial y}{\partial \eta}\right)^2\\ b=\left(\frac{\partial x}{\partial \xi}\right)\left(\frac{\partial x}{\partial \eta}\right)+\left(\frac{\partial y}{\partial \xi}\right)\left(\frac{\partial y}{\partial \eta}\right)\\ c=\left(\frac{\partial x}{\partial \xi}\right)^2+\left(\frac{\partial y}{\partial \eta}\right)^2 $$

The form of the governing equation is the same as we solved in the elliptic grid-generation method! This means we can use the functions we wrote previously for that purpose.

Boundary Conditions in the Computational Domain

Though the form of the governing equation is the same as the PDEs for the elliptic grid-generation method, the boundary conditions are not. We now examine the boundary conditions.

On the outer boundary, we expect the flow velocity to be equal to the free-stream velocity, that is:

$$ \left\{ \begin{array}{l} u=\frac{\partial \Psi}{\partial y}=V_\infty\cos(AOA)\\ v=-\frac{\partial \Psi}{\partial x}=V_\infty\sin(AOA) \end{array}\right.,\ on\ outer\ boundary $$

where $AOA$ represents angle of attact, and $V_\infty$ is the free stream velocity.

Since $V_\infty$ and $AOA$ are constant, we can simply integrate the above boundary conditions for the stream function on the outer boundary

$$ \Psi=-xV_\infty\sin(AOA)+yV_\infty\cos(AOA)+C_1,\ on\ outer\ boundary $$

In potential flow, the values of $\Psi$ are not important; what matters are the derivatives of the stream function. Since the constant $C_1$ in the above equation will disappear in the derivatives, we can simply make $C_1=0$ here. Now we can apply Dirichlet boundary condition on the outer boundary:

$$ \Psi_{i,j=-1}=-x_{i,j=-1}V_\infty\sin(AOA)+y_{i,j=-1}V_\infty\cos(AOA) $$

Next, we examine the boundary conditions for nodes on the airfoil. For walls in potential flows, the no-slip boundary condition does not apply, given the lack of viscosity. The boundary condition for walls in potential flows is that the flow can not penetrate the wall (no-thru BC). This implies that the wall itself is a streamline.

$$ \Psi=C_2,\ on\ the\ airfoil $$

or

$$ \Psi_{i,j=0}=C_2 $$

Since we have determined the values of $\Psi$ on the outer boundary, we cannot determine the value of $\Psi$ on the airfoil, $C_2$ (because we don't know which streamline on the outer boundary corresponds to the streamline on the airfoil surface). Fortunately, we will resolve this problem later after we introduce the Kutta condition.

Finally, the boundary condition on the dividing line $\overline{AC}$ is again given from periodic boundary condition.

$$ \Psi_{i=0,j}=\Psi_{i=-1,j} $$

The Kutta Condition

In potential flow, due to the lack of viscosity, the flow under an airfoil will flow back to the upper surface after it has passed the trailing edge, even at a low angle of attact and low free stream velocity. This is not physical. The Kutta condition ensures this situation won't happen: it makes the flow smoothly pass by the trailing edge. More details can be found in any textbook covering inviscid flow or potential flow. Here we only introduce the implementation of the Kutta condition in our code.

According to reference [3], for the airfoils with finite-angle trailing edge, like the NACA 2412, the Kutta condition results in the value of $\Psi$ on the node next to the trailing edge being the same as that on the trailing edge, in order to provide a smooth streamline near the trailing edge. In other words,

$$ \Psi_{i=0, j=1}=\Psi_{i=0, j=0} $$

Recall the boundary condition on the airfoil surface: $\Psi_{i,j=0}=C_2$. We can determine $C_2$ and $\Psi_{i,j=0}$ now. In each iteration when solving the governing equation, let

$$ \Psi_{i,j=0}=C_2=\Psi_{i=0,j=1} $$

The problem of the boundary condition on the airfoil surface is then resolved. We can start to solve for the stream functions!

Solve It Now!!!

Let the angle of attack be $15^\circ$ and free stream velocity be $70\ m/s$. Other parameters are the same as we used in previous sections.

Though the values of stream function are not important, conventionally, we make the streamline on the airfoil surface equal to zero.

Finally, the streamlines are shown in the following figure.

We successfully solved the potential flow over a NACA airfoil. We can start to try other parameters now!

Reference

[1] NACA airfoil, Wikipedia, http://en.wikipedia.org/wiki/NACA_airfoil

[2] K.A. Hoffmann, and S.T. Chiang, Computational fluid dynamics, Vol. 1, Wichita, KS: Engineering Education System (2000).

[3] F. Mohebbi and M. Sellier. On the Kutta Condition in Potential Flow over Airfoil, Journal of Aerodynamics (2014).